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Abstract 

Some ideas studied by Lele ( 1993ll . under a Gaussian perturbation model, are gen¬ 
eralised in the setting of matrix multivariate elliptical distributions. In particu¬ 
lar, several inaccuracies in the published statistical perturbation model are revised. 
In addition, a number of aspects about identifiability and estimability are also 
considered. Instead of using the Euclidean distance matrix for proposing consis¬ 
tent estimates, this paper determines exact formulae for the moments of matrix 
B = , where X*^ is the centered landmarks matrix. Consistent estima¬ 

tion of mean form difference under elliptical laws is also studied. Finally, the main 
results of the paper and some methodologies for selecting models and hypothesis 
testing are applied to a real landmark data, comparing correlation shape structure 
is proposed and applied in handwritten differentiation. 


1 Introduction 


Statistical theory of shape has emerged as one of the most versatile techniques of classifi¬ 
cation and comparison of “objects” in a number of disciplines. By its theoretical nature, 
the matrix multivariate distribution analysis fits very well into the shape analysis, but at 
the same time have involved strong open problems on estimation of location and scale pop¬ 
ulation parameters based on the exact distributions, forcing the application of several less 
robust approaches, which were considered appropriate at first, but lat er received important 
critics from different experts on morphometries and related fields, see Eeii (Il993h . 


Among the addressed lacks we can cite the use of asymptotic distributions, tangent plane 
inference, isotr opic models, Gaussian assu mptions, and procrustes theory; for details of such 
techniques see iDrvden and Mardial (|l998h and the references therein. 
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Now, some attempts have been published recently avoiding the above restrictions and 
considering inference via likelihood function using the exact shape distributions, the new 
theory was termed generalised shape theory by hnding the exact shape densities indexed 
by families of distributions of elliptical contours. According to the geometrical filters on 
shape, the resulting exact invariant distributions are expressed in terms of series of func¬ 
tions termed Jack polynomials, which w ere uncomputable for decades, and only recently with 
the works of iKoev and Edelman ( 2006ll . the individual polynomials could be computed but 
series of them have involved serious problems in inference of population parameters via like¬ 
lihood method. A number of approaches with a meaningful computational success in the 


position, see Goodall and Mardia (Il993h and Dfaz-Garefa and Garo-Loperal (20141. sinmi- 

lar value decompositions, see Le and Kendalll ( 

1993h. Goodall (19911. Dfaz-Garefa et a/. 

(20031. Dfaz-Garefa and Caro-Looera ( 

2012al and Dfaz-Garefa and Caro-Looera (2012tJll. 

affine see Goodall and Mardia 1 ( 1993ll. 

Dfaz-Garefa et al. (2003^. Caro-Looera et al. (2009ll. 


Dfaz-Garcia and Caro-Lopera " ([2013 1. However, a feasible approach dealing with com- 


see 

putable exact densities and a likelihood function based on polynomials of very low degree was 
published recently, letting robust estimation on location and scale population parameters 
very accurate; it models shapes under certain conditions via affine transformation, which 
means that it removes from objects, any geometrical information of rotation, translation, 
scaling and uniform shear. Meanwhile, the similarity (Euclidean) transformations via QR, 
SVD, Pseudo-Wishart (invariant under rotation, translation, scaling) capture the attention 
of most of the users of shape theory and is the source of the main critics. 

Under Euclidean transformations the shape distributions are extremely difficult to com¬ 
pute and then the associated inference, it forces the use of isotropic models, an assumption 
which is unrealistic in biology, for example, since it says that landmarks vary indepen¬ 
dently of each other along different axes but are correlated along a fixed axis. In fact, 
this isotropy assumption is very common in literature, leaving the problem of testing solely 
whether shapes are equal; but for biologists, for example, they want to identify the cor¬ 
relation structure of landmarks and the structures of shape which absorbs the meaningful 
differences. Moreover, estimation of a full covariance structure would be the desirable result, 
because, correlation among landmarks and axis is important, but the correlation among ob¬ 
ject s^n__the_sam£ffi_shouldprwide a complete comprehension of the involved populations, 
see 


Lele and Richtsmeier (199^ 


Instead of estimation via likelihood method, some authors have proposed, the Gaus- 
the method-of-moments estimators of the mean form and the variance-covariance 


Sian case 


struct u re wh i ch are consistent and simple to c ompute, see for example iLele and Richtsmeier 
(1991), Leld ( 1993 1. Richtsmeier et all 1 2002ll and the references therein. In fact, the tech¬ 
nique was set as a critics of generalised procrustes analysis, by proving that application of 
the last analysis yield inconsistent estimators of the mean form, mean shape, and variance- 
covarian ce struc t ure, a nd then all the statistical inference proc edure s can p roduce inaccurate 
results. IWalkerl ( 200 ill recently reiterated the conclusions of Leld ( 199311 by reporting the 
inability of Procrustes methods to estimate the correct variance-covariance structure and 
the associated implications for statistical inference. This aspect, is crucial because Pro¬ 
crustes analysis is one of the most common method of estimation in several fields such as 
morphometries. 

Given the computation problems of maximum likelihood estimators, the method-of- 
moments estimators emerges as one of the promissory techniques in shape theory, however 
some open p roblems mu st be studied deeply. For example, some inaccuracies of this model 
presented in Leld (1993), assuming a matrix multivariate Gaussian distribution must be 
nuanced first, and second, the method should allow non Gaussian samples, a realistic and 
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very common problem in morphometries and the usual applied areas for shape, a suitable 
solution comes from families of elliptical contoured distributions, which exhibit lighter or 
heavier tails, or greater or less kurtosis than the Gaussian model. Setting generalised shape 
theory also must include criteria for selecting models and hypothesis testing, in order to 
provide an integrating theory suit able t o be a pplied in meaningful scenarios. 

Clarifying the inaccuracies of Le^J 19^ ) id eas, and their con necti on with some theo- 
retical studies bv iMagnus and Neudecker ( 1979ll . Muirheadi ( 1982h and Diaz-Garcial (1994) 
should give a unihed theory setting the isolated Gaussian approach into the general frame¬ 
work of the existing generalised matrix multivariate elliptical shape theory. 

Thus, estimation of mean form and mean form difference under elliptical laws is placed 
in this work as follows: Section [5] clarifies some results of the published Gaussian case and 
propose the generalisation in the context of matrix multivariate elliptical distributions, it 
includes some properties of matrix multivariate elliptical distribution, identifiability and 
estimability of the parameters of interest, the perturbation model under a matrix multivari¬ 
ate elliptical distribution, and invariance and nuisance parameters. Then Section [3] studies 
the consistent estimation of the population parameters under dependence and independence 
and provides exact formulae for the moments estimators. Section 0] provides a consistent 
estimation for a general non-negative definite correlation matrix. The analysis also includes 
extensions to elliptical models of form difference under the perspective of Euclidean Distance 
Matrix, see Section [5] Finally, a complete example collecting the main results of the paper 
and proposing some selecting model criteria, is proposed in Section [51 


2 Preliminary results 


In this section we review some notation and distributional results. Also, the statistical 
model to be used throughout the pa per, is esta blished and analysed. In particular some 
inaccuracies of this model presented in lLelel (19931), assuming a matrix multivariate Gaussian 
distribution are corrected and then is generalised to the case where an matrix multivariate 
elliptical distribution is assumed. 


2.1 Matrix multivariate elliptical distribution 

A detailed discussion of the matr ix m ultivariate elliptic a l dist ribution can be found for 
example in Fang and Zhang! ( 199d) and Gupta and Vargal (1993), among many others. 


Remark 2.1. For matrix multivariate Gaussian and elliptical distributions, traditionally 
are used two forms for establish that a random matrix Y has one of these distribution. For 
example in matrix multivariate Gaussian case, this fact is written as 

Y^A/'iCxD(M,S,©), 


see 


Arnold ( 198lll . lDutilleu 



However, as is study in lLeL 
S and © are not identifiable one-by-one, but 


and Fan g and Zhand ( 199Cl) among many others authors. 
iDutilleul (1999) among others, in general the parameters 
S 0 © or © 0 S are identihable. Here 0 
denotes the usual Kronecker product. In addition, given that Gov(vec Y) = © 0 S, and 
Cov(vecY'^) = S0©, many other authors use the notation 


Y ~ A/’/cxd(A‘, 51 0 ©), 


where “vec” denotes the vectorisation operator, see lMuirhead ( 19821) and Gupta and Varga 
(|l99311 . Analogous situation is present for matrix multivariate elliptical distributions. We 
shall use this last notation. 
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Definition 2.1. It is say that Y has a matrix multivariate elliptical distribution, with 
location parameter matrix fi G and scala parameter matrix S ® 0 G 

where S is a definite positive matrix, S > 0 and © > 0, with S G and © G 

if its density function with respect to Lebesgue measure is given by 

dF,{Y) = \'E\-^/^\&\-^^^h[tT&-\Y - - ti)]{dY), (1) 

where the function ft. : 5i —?► [0, oo) is such that u^^^^~^h{u)du < oo. The function ft is 
termed the density generator. Its characteristic function is given by 

V’y(T) = etr(i/x'^T)^(trT©T'^S), (2) 


with i = \/^, (j) : [0,oo) —>■ 3? and etr(-) = exp(tr(-)). This fact is denoted as Y ^ 
SKxoi/J-, S ® ©, ft). In addition, observe that the characteristic function exist still when S 
and/or © are semidefinite positive matrices; in such case it say that Y has a singular matrix 
multivariate elliptical distribution, see Remark Ti. 21 below. 

In addition, note that Cov(vecY) = cq© 0 S, and Cov(vecY^) = cqS 0 © where 

CO = -2</'(0), 


</'(«) 


d(j){t^) 

dt 


t=o 


see ( Fang and Zhane . I990l Theorm 2.6.5, p. 62) and ( Guuta and Varga . I993L Corollary 
3.2.1.1, p. 94 and Theorem 2.4.1, p. 33). 


Is easy to see that if 



( (Til 0-12 

(721 <722 

0\K \ 
<^2K 

and © = 

/ 011 012 

021 022 

01D \ 
0211 

\ Cifi OKI • 

<7KK / 


V 0r>i 0r>2 • 

0fcn / 


Then from (|Fang and Zhangl . Il990ll 


1- ^(i) ; <7ii©, ft), i — 1 , 2 ,..., K, 

2. Yj £K{fJ.j,0jj'E,h), j = l,2,...,D, 


this is. 


1. Cov(Y(i)) = coCTi*©, i = l,2,...,K, 

2. Cov(Y,) =co0jyS, / = 1,2,...,D. 


These two last affirmations are incorrect stated in iLelel (| 19931 1 in the context of the pertur¬ 
bation model. For this asseveration observe that this class of matrix multivariate elliptical 
distributions includes Gaussian, contaminated Gaussian, Pearson type II and VII, Kotz, 
Jensen-Logistic, power exponential and Bessel distributions, among others; these distribu¬ 
tions have tails that are more or less weighted, and/or present a greater or smaller degree 
of kurtosis than the Gaussian distribution. 
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2.2 Identifiability and estimability of the parameters of interest 

Now some aspects about the identifiability and estimability of the parameters (/i, S ® 0) 
are studied. 

Note that the density © can be write as, see dMuirhead 1 . 119821 p. 79) and ( iGuota and Vargcl 


19931 Theorem 2.1.1, p. 20), 


dF ^{vecY'^) = |S(8)©r^/^/i[vec'^(Y-/x)(S(g)©)-%ec(Y-/i)](dvecY^), 


(3) 


using the fact that vec^ X(DB 0 C^)vecX = tr(BX^CXD), with ve c'^X = (vec X)^ , 
and that for A € an d B g then |Ar |B|" = |A 0 B|, see ()Muirhead I. ll982L 

Section 2.2, pp. 72-76) and ( Fang and Zhand . 1990l Section 1.4, pp. 11-13). Then, denoting 
vec Y^ = y G and S 0 0 H, the density ([3]) define the distribution of the vector y; 

moreover, y ^ £lif£)(vec/x, H, h). 

Now, assume that our data consist of a sample of matrices of size n from a given popu¬ 
lation, namely Yi, Y 2 ,..., Y„, and define the random matrix 


Y = (vec Y;), vec Y 2 ,..., vec Y„) £ 


ynxKD 


From iDiaz-Garcial (Il994ll . assuming that Yi, Y 2 ,..., Y„ be independent, the density func¬ 
tion of Y admit the expression 

dF^(Y) = |S 0 0|-"/2/i[tr(Y - M)(S 0 0)-i(Y - M)'^](dY), (4) 

where 

M= l„vec^/x G 

and 1„, = (1,1,..., 1)^ G this is, Y ^ S 0 0 0 I„, h). Thus, taking p = KD 

in ( Fang and Zhanj . Il990l Theorem 4.1.1, p.l29), and given that KD < n and h{-) being 
nonincreasing and continuous, we have that the maximum likelihood estimate of (vec /x, S 0 
0 ) is 

(vecfi, S 0 0^ = (y, AinaxS), 

where Amax is the critical point where the function h*(X) has its maximum, with 

h*{X) = X-^^'^/^h{KD/X), 


1 , 


y = -Y" 1„ G and S = Y" H„Y G » 


KDxKD 


1 


where H„ = I„-Inll, defines an orthogonal projection, this is, H„ = H'E = Hi . Or 

n 

alternatively 

.. n n 

y = - X! Yf, and S = ^(vec Yf - y)(vec Yf - y)^. 


from where the estimator of is 


M = Y=-^Y.. 

71 


2=1 


From ( Fang and Zhanj . 1990l Section 4.3), several properties of the maximum likelihood 
estimators J1 and S 0 0 are obtained as: sufficiency, completeness, consistency and un¬ 
biasedness. Specifically, for Y with the finite 2nd moment and h{-) be nonincreasing and 
continuous, 

p. = Y and S 0 0 = 


and S 

are unbiased estimators of /x and S 0 0. 


2(l-n)^/>'(0) 
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Remark 2.2. Observe that when Y ~ £kxd{iJ‘, S ® 0, /i) and its columns and/or its rows 
are dependent linearly, is say that Y has a singular matrix multivariate elliptical distribution. 
Then Y has density with respect to Hasusdorjf measure. Moreover, such dependent linearly 
among its columns or its rows is archived in the rank of S and/or © matrices and is 
denoted as: Y ^ n ( fx, S Q &,h), where s = ra nk(S) < K and r = rank(0) < 
D, see ( Gupta a nd Varga, Il993 , Definition 2.1.1, p . IQl lPiaz-Garcia and Gonzalez-Farias 


( 2005 1 and Diaz-Garcia and Gutierrez-Jaimej ( 2006ll . As in the singular matrix multivariate 
Gaussian case, the maxi mum lik e lihoo d esti mators in si ngular matrix multivariate elliptical 
models remain valid, see Khatri ( 1968 1 and ( Rad . 19731 Section 8a.5, pp.528-532). 


2.3 Perturbation model under a matrix multivariate elliptical dis¬ 
tribution 

LetX G a random matrix representing the geometrical figure comprising K landmark, 

or labeled, points of dimension D , such that K > D. This matrix X is termed landmark 
coordinate matrix, see II^ (I 1993 I 1 . 

Let Xi, X 2 ,..., X„ be a independently sample of size n of landmark coordinate matrices 
X, G i = 1,2,... , n, from a given population. 

The statist ical model t o be considered in this work is a generalisation of the perturbation 
model used bv iLele ( 1993 1 among others authors. Let G corresponding to the mean 

form. Let 

Xi = (/XEdTi-|-ti, , z = 1, 2,..., n, (5) 

where Ei ^ £kx are orthogonal matrices representing rotation 

and/or reflection of (/i-f E^), and G are matrices such t hat = Ifcaf representing 

translation, for some G From ( Fang and Zhanj . ll99Cll eq. (3.3.10), p. 103) or 

( Gupta and VargJ.ll993l Theorem 2.1.2, p. 20) we have 


Xi ~ + ti, Sif 0 tJ z = 1,2, 


,n. 


( 6 ) 


Parameters of interest are {pL,T,K ® ^d) and (r/’,ti) z = 1,2, ...,n a re th e nuisa nce pa¬ 
rameters. An detail explained of this perturbation model is given in iLelel (|1993[1 among 
others. 

Alternatively, the model ([S]) can be write as: 


with 


vecX^ = diag(G)vec( 

/ Ix 


- E)^ -(- vecT^ 




0 


0 


l-K I 


rl’ 


diag(G) = 

0 0 

or the model ([S]) may be rewritten in the form 

X = diag(M -f E)G'^ -H T, 

with 


LK ' 


\ 

) 


diag(M -I- E) = 


( vec"^)/!-1-Ei)^ 0 

0 vec^(/i-I-E 2 )^ 


0 

0 


vec^(/2-I-En)"^ / 
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where 



/ vec^ Xf \ 


/ vec^ Ef \ 


/ vec^ t|’ \ 


vec^ 


vec^ El’ 


vec^ t|’ 

X = 


, M = 1„ vec^ /j,^, E = 


T = 



V vec^Xj ) 


1 vec^E^ y 


V vec^ti; y 


and G = (Ik ® ® r|’| • • • \Ik ® T^) where 

E ^ £nxKD{ 0 ^ In ® 0 '^D, h), 

or 

vecE^ ~ f„K£)(vecO,I„ ® ® S^i, /i), 

Hence 


vecX^ ~ fn/CD (diag(G) vecM^ + vecT^,diag(G)(I„ 0 0 S^i) diag(G)^, h) 

Note that, recalling that for x and y vectors, vec(yx’^) = x 0 y, then 


diag(G) vecl 


lif 0 0 

0 lif 0 rl’ 


\ 0 0 

^ (Ik 0 rf’) vec/x^ ^ 
(Ik 0 r|’) vec /x^ 


V (Ik 0r^)vec/x^ / 
/ vec(/xri)'^ \ 
vec(/xr2)'^ 


\ vec(/xr„)'^ y 

and diag(G)(I„ 0 Sk 0 Sk) diag(G)^ is 


y Sk 0 rf’SKTi 

0 Sk ' 


- 2 2 


= SeJ 


Sk 0 rf SkFi 


i=l 


0 \ 
0 


(l„ 0 vec /x^) 


•■if ' 


y 


0 

0 


SK0r^SKrji y 


where if e" the ith column unit vector of order n, then E"j = e'*(e")^. 
Finally observe that 


E(vecX^) 


/ vec(/xri)^ \ 
vec(/xr2)^ 


+ vecT^ 


\ vec(/xr„)^ y 


^ ey 0 (vec(/xrj)^ + vectf) 
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Then 


EOQ {vec{fj,Tif + vec tf)' 


Therefore 


^nxKD 


Y, er (vec(/xr,)^ + vec tf)"’, ^ E" 0 ® tJ S^T,, h . 


2.4 Invariance and nuisance parameters 

In general, when a model contains nuisance parameters, the hrst step i s to remove them. 
As in the matrix multivariate Gaussian model considered bv iLelel ( 1993li . under an matrix 
multivariate elliptical model this objective is achieved through a simple transformation. 
From dH) 

Xi ~ EKxoifJ-Ti + ti, S/f 0 rf h), i = 1,2,... ,n. 


Recalling that = 0^ and = 0^, then, defining X° = we have 


X? 


'-K 

SJ, 0 Tj ^dT., h), i = l,2,...,n. 


(7) 


where fi* = Hif/x and = H/f , Hi^t^ = H/f l^af = 0 for alH = 1, 2 ,..., n, and 

/.i* is such that its columns sum to zero, that is, it is a c entered matrix . _ 


Given t hat K > D and that rankfSt^j = A T — 1, f rom iDiaz-Garcia and Gonzalez-Farias 


( 2005 1 and Diaz-Garcia and Gutierrez- Jaime3 (l2006ll we have that 

= ^grwl{D,^*^,'ED,n,h), ^ = l,2,...,n. (8) 

where 

q = min((Ar — 1), D) and A“ is any symmetric generalised inverse of A such that AA A = 
A = A"^. This is, has a generalised singular pseudo-Wishart distribution, which is 
independent of noise parameters. 

Remark 2.3. Observe that can be write as 

B, = X=(rfS^R)-i(X?)^ = X?rfS-iR(X=)^ = 

where Y^ = X^F^ and is such that 

Y, ^K^^D,h), i = 1, 2,..., n. 

In particular if and 

x= = (X5 JXL|...|x^o 


with 


XS 


we have that. 


furthermore, 


d=l,2,...,D-, i = l,2,...,n, 

D 


d^l 


B, - gVW],{D, X*J^, Id, n,h), 1 = 1,2,..., n, 

where 


(9) 














Remark 2.4. The result in iLelel ( 1993[l is obtained a s particular case of ([5]), with the 
difference that the matrix of noncentrality paramet er in ( 1993^ is defined as /x*(/x*)^ 
and we use Jl = notation used in ( Muirhead Il982l Definition 10.3.1, pp. 

441-442). 

In addition, defining X'^ as X we have 

vec (X")^ = [I„ 0 (Hfe 0 Ip)] vecX^ 

hence, X'^ = X(Hfc 0 Ip). Now, observing that (H^ 0 Ip) vectf = 0, for alH = 1,2,... ,n. 
Then 


n,(K-l)D 

nxKD 


where = HfeS^-IIfc. 


fcSR- Jlfc 


e- vec^(M*R)^, ^ 0 0 rf SpR, h 

i=l / 

As in held ( IQQSll . assuming that Sp = Ip, and recalling that 

n 


2=1 


we have 


^n,{K-l)D 

■'nxKD 


er vec^(M*R)^, I„ 0 0 Ip, /i . 


( 10 ) 


3 Consistent estimation of fx and 


Alternatively to the use of the Euclidean distance matrix showed in Lele ( 1993h with the 
aim to propose consistent estimations, we use directly the first two moments of the matrix 
B with the same object. 

When is considered a model where the perturbation of landmarks along the D axes are 
independent and identical to each other, formally we are assume that Sp = Ip under a 
matrix multivariate Gaussian case. However, this same assumption is not to hold in matrix 
multivariate elliptical case. Under a matrix multivariate elliptical case is possible to consider 
two cases: 

1. Independence and not correlation among landmarks and 

2. Probabilistic dependence and not correlation among landmarks. 

In both cases Sp = Ip and the moments of matrix B are different in each case. 

Remark 3.1. Recall that under matrix multivariate elliptical distribution, only in the 
Gaussian case the not correlation and independence are equivalent. Then suppose that the 
vector Z = (zi, 22 )^ has a bi-dimensional elliptical distribution and Cov(Z) = I 2 then Zi and 
Z 2 are independent if and only if Z has a bi-dimensional Gaussian distribution. But if Zi, have 
a uni-dimensional elliptical distribution for i = 1,2, and Var( 2 :j) = 1 a nd Covf^i, Z 2 ) = 0 , 
Zj, i = 1,2, are not correlat ed and can be considered independent, see (|Gupta and Varga . 
I 993 L Section 6.2, p. 1) and ( Fang. Kotz and n 3 . 1990L Section 4.3, p. 105). 


Summarising, given 


B = YY^=^ydyJ, 


( 11 ) 


next, we find the first two moments of B assuming that Sp = Ip, i.e. when y^: a) are not 
correlated and independent; and b) are not correlated and dependent. 
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3.1 Moments of B under dependence 

By completeness initially we assume that ^ and for convenience denote = ©, 
= S and /i* = /x. 

With this goal in main, suppose that Y ^ {fi, S ® 0, h), with 

Y = (yi|y2|-" |y£>) and /x = (/Xi|/X2|-- - l/x^,). 


Observing that for x, y € 5ft", vecxy^ : 
vecyy"^ vec^ yy"^ = y O y"^ O y ® y^, see 


= y (H) X, xy^ = x 0 y" = y" 
Magnus and Neudecker ( 1979h . 


T 

y = 


O X and thus. 


Theorem 3.1. Let Y ^ £’j^^^’^(/x, S ® 0, h). T/ien 


1. E{YecY vec^ Y) = co(0 ® S) + vec/xvec^/x, 
and ElvecY vec Y^ ® vec Yvec Y^) is 

= ko[(I(xd )2 + Kkd)(© ® S (g) © (gi S) + vec(© 0 S) vec^(© 0 E)] 
+ Co (Ij^2 + Kk ) [vec /X vec^ /x 0 (© 0 S) + (© 0 S) 0 vec /X vec^ /x] 
+ co[vec(© 0 E)(vec^ /x/x^) + (vec^ /x/x^) vec(© 0 S)] 

+ vec/x vec^ /X 0 vec/x vec^ /X, 


where is t/ie com mutation matrix, see Magnus and Neudeckei 1 197^) . and cq = E{u^) 

and 3ko = E{u'^), see iGuvta and Varaa . 199 A . p. 127), 


E(U-) = i 


dt^ 


a«<i £(„*)= 


t=0 


dt^ 


t=0 


Where '0u(t) = (/>(t^) is the characteristic function of univariate elliptical distribution. Some 
particular values of co and kq, are summarised on Tabled 


Proof. This is obtained differentiating ([2]) and observing that, see lDiaz-Garcia and Gutirrez Jimez 


E{\ecY i^ivec^Y) = 


and 

E{vecY 0 vecY^ 0 vecY 0 vecY"^) = EfvecY vecY^ 0 vecYvec Y^) 

_ j__ 5'^^vecY(vecT) _ 

i"^ (3 vec T9 vec T^ 5 vec T9 vec T^ vecT=o' 

□ 


Efvec Y vec^ Y) 

1 92V’vecY(vecT) 


i^ 9 vec T9 vec T^ 


vecT=0 


Now, given 


we have 


B = YY^ 




£;(B) = E (YY^) = E 



D 


Y. ^ (y^Yd) ■ 
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Table 1: Particular values of cq and kq. 


Distribution 

Co 

Ko 

Multiuniforme_2 

1 

1 

3 

Gaussian_^ 

1 

1 

Kotz^ 

p 2J\I -\-l 

^ L 2s J 

L 2s J 

1 /,,-p r2Ar-i] 

' ^ . 2s . 

3r2Ar 


ra 


m — 2 

(m — 2)(m — 4) 

Pearson Type IF 

1 

1 

2m + 3 

(2TO + 3)(2m + 5) 

Pearson type Vllf 

m 


to 

1 

CO 

(2N - S)(2N - 5) 


(Fang. Kotz and N dll990l . Theorem 3.3, p. 72). 
i cSta and Vargal . [l993l . Remark 3.2.2, p. 125). 

( Nadaraja^HooB)' ) where r, s > 0 and 2N + 1 > 2. 

Gupta and Var^. |l993|. p. 128), or l|Fang. Kotz and Nd . [ToSol . p. 88), where m > 0. 


“From 
*'From 
“From 
"^From 
“From 

^From I Fang. Kotz and Nd . 


I ^u£t^^n^^arga|Jiyy^j_p. 

I Fang. Kotz and Ngl 199ol. ! 


1991 


Section 3.4.2, p. 
Section 3.3.4, p. 


89), where m > —1. 

84), where N > 1/2, m > 0. 


And remembering that for Y G in general 

Cov(vec Y) = E{vecY vec”^ Y) — E{vecY)E{vec^ Y). 

Therefore 

Cov(vec B) = Gov (vec (YY^)) = Cov l^^vec (ydyj)^ 

vec {ydyj )) vec^ (y^y' 

\d=l ) \s=l 

“ ^ (y^yl)^ ^ (y^yJ) 

1 

^ E (y^yj) 


= E 


.d=l s=l 


\d^l 


— vec 


vec 


^ E (y^yf) 




vs=l 


( 12 ) 


Then, we need to find E (ydyj) and E (y^yl 0 y^yl) ■ These moments are obtained in the 
following result. 

Theorem 3.2. Assume that Y ~ S 0 0, h), with 

Y= (yi|y2|---|yD) andn = (Ati|/X2|---|MD), 

and 0 = (Ods)- Then 
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1- E{ydyJ) = + Hdfj-J. 

2. And 


E {ydyj'^yd.y’^) = ko^'L[(Ik=+K; if)(S0S)+veci:vec'^S] 

+ coOds[{lK^ + (g) S + S 0 Md/xf)] 

+ coOds [vec S vec^ fJ-dt^I + vec HdfJ-I vec"^ S] 


Proof. The results is obtained as consequence of Theorem 13.11 observing that yd = Ye^, 
then 

EiydyJ) = E{vecydvec^ yd) = E{vecYe° vec^Ye^) 

= (ef ^ 0 Idf)ii;(vec Yvec^ Y)(ef 0 lif) 

= (ef ^ 0 ldf)(co(0 0 S) + vec/x vec^/x)(ef 01^) 

= CoOdd'^ + fJ-dP-d ■ 

This least result is obtained noting that, vec ABC = (C^ 0 B) vecB, a 0 A = oA and 
(A 0 D)(B 0 E)(C 0 F) (ABC 0 DEF). Similarly, 

E (y^yj ® ydyj) = R^E(vec Y vec Y”^ 0 vec Y vec Y^)Ri 

with R^ = (e^ ^ 0 lif) 0 (e^ ^ 0 lic) and Ri = (ef’ 0 Ik) 0 (ef* 0 Ik)- The desired result 
is obtained observing that: for A g and B € Kmn(A 0 B) = (B 0 A)Kts and 


that Kmm = Km see lMagnus and Neudeckei ( 1979f) . 
Consider the following definition. 

Definition 3.1. Let A £ such that 



f An 

Ai2 

-^In 

_ 

A 21 

A 22 

^2n 


^ A^i 

Am2 

A 

■^mm 


then 




A„- £ 5R’' 


□ 


771,n 

ffl 


i=l d=l 


m,m m 

If m = n then, ffl =ffl. 

ij ij 

In addition, let A = ( An) and B = (B^-) partitioned matrices. Then if 0 denotes the 
Khatri-Rao product, see ( Rad . 1973 . p.30), 

A 0 B = (Aij 0 Bij)^^ . 

In particular, note that if C = (cy), then 


Moreover, 


C 0 A — (cijAij)^j . 


(C0A)=EE("bA,d), 

^ j 
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Theorem 3.3. Suppose that Y ~ S (g) 0, h), with 

Y = (yi|y2| • • • lyn) and /r = ■ ■ ■ |Mu). 


And define 


Then 


And 


B = YY^ 


'^ydyj- 


i?(B) = Co tr(0)S + . 


Cov(vecB)) 


(I^2+KK){Kotr(02)(S^ 


Co 


D 


1,3 

D 


(0 0 vec fi vec^ /r) 0 S 


+S0 ffl (0 0 vec fj, vec^ /x) 


+ [^otr (0^) — Cotr^(0)] vecS vec^ S 


Co I vec S vec^ ffl (0 0 vec /x vec^ /x) 


D 


vec 




(0 0 vec /X vec^ /x) vec^ S 


+ tr(0) [vec S vec'^ /x/x^ + vec vec S] } . 
Proof. This is a consequence of dD, (ini), Definition [XT] and Theorem 2. 
Corollary 3.1. In Theorem \3.‘A assume that ® = Id. Then 

E(B) = Dcq'E + /x/x^. 


And 


□ 


Cov(vecB)) = (1^2 + Kif) {Dko(S 0 S) + Co [/x/x"^ 0 S + S 0/x/x"^] } 

+ D [kq — -Dcq] vec S vec^ S 
+ (1 — D) Co [vec S vec"^ fifi^ + vec /x/x"^ vec^ S], 


In un ivariate case, when p, = 0, these results were obtained in general and a particular 
cases in ( Quota and Varg^ll993l Theorem 3.2.13 and Example 3.2.1), with a several minor 
errors. In particular, for general case they write D^{kq — Cq) and for matrix multivariate T 
distribution they write D{ko — Cg), with D = n — 1, instead of D (kq — Dcq). 


3.2 Moments of B under independence 

Let Y and fi such that 


Y = (yi|y2|-- - |y£>) and /x = (/X 1 I/X 2 I • • • l/x^,), 
where yi, y 2 ,..., y_D are independent and 
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and by independence, Cov{yd, y^) = 0 , for d ^ s = 1, 2 ..., D. 
Given 

D 


B = YY^ = ^ 


ydyj, 




we have 


i?(B) = (YY^) = i? y.yj] = i? (y,yj) . 

\d^l / d=l 

And under assumption that d = 1, 2,..., Z) are independent, 

Cov(vec B) = Gov (vec (YY^)) = Gov (e vec (ydyj)^ 

D D 

= ^Gov(vec(ydyJ)) =^Cov(yd(g)yd). 


d=l 




Then, we need to find E (y^yj) and 

Gov (yd (g) yd) = E{{ydYd) {VdYd)'^) - E {ya 0 yd) E {yd Yd)'^ 

= E (ydyj (g) YdyJ) - E (vecydyj) E (vec^ ydyj) ■ 

These results are obtained in the following 

Corollary 3.2. Let yd £^j^~^\fJ.d^^dd'^;h), d= 1,2,. where yi,y 2 ,... ,yD are 
independent. Then 

1- E{ydyJ) = co6»ddS + Md/^J- 

2. And Cov(yd g) yd) = Gov(vecydyJ) 

= (Idf 2 + K.k) {Koddd(^ g S) + co9dd [MdMd g S + S g /XdMd] } 

+ ^ddi^o — Cq) vec S vec^ S. 


Proof. It is follows from Theorem 13.21 taking d = s. 

Theorem 3.4. Suppose that yd ^ £^~^\p.dy^dd'^',h), d = 1,2,..., D, with 
Y = (yi|y 2 | • • • lyu) and n = (/X 1 I/I 2 I • • • Imu), 


and let © = diag(dii, 022 , ■ • ■ j 9dd), and 


D 


B = YY^ = 


ydYd- 


Then, 


E{'B) = Co tr(©)S + 

Gov(vecB) = (Id ^2 + Kdc) {kq tr (®^) (51 g S) 

E ^ddP'dP'dj g S + S i 


Co 


i=l 

2 


+ {kq — Cq) tr (©^) vec S vec^ S 


D 


E ^ddtJ-dP-d 


\d=l 


□ 
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Proof. From Corollary 13.21 


Similarly, 


E{B) = E{YY^) = E(Y.ydyd] =J2^iydyd) 

\d^l / d^l 

D 

= i^oOdd's + fj-dfdi) 

d^l 

D 

= cotr(0)S + = cotr(0)S + 


Cov(vecB) 


Cov (vec (YY^)) = Cov 
^Cov(vec (ydyl)) . 


D 


vec (ydyj) 




from the desired result is obtained. 


Now if 0 = I^), we have the following results. 

Corollary 3.3. Let yd (/x^, S;/i), d = 1,2 ,where yi,y 2 , 

independent. Then 

1- E{ydyJ) = CoS + /Xrf/xJ. 

2. and Cov(yd 0 yd) = Cov(vecydyJ) 

= (Iif2 + Kx)[ko(S 0 S) + co(/x^/xJ 0 S + S 0 /x^/xj)] 
+ (ko — Cq) vec S vec^ S, 


Proof. It is immediately. 

Theorem 3.5. Suppose that yd ^^(/x^, S;/i), d = 1,2,... ,D, with 
Y = (yi|y 2 | • • • lyn) and /x = (/X 1 I/X 2 I • • • l/x^,), 


and let 


Then, 


E{B) 
Cov (vec B) 


B = YY^ 


Yyd-y'd- 


HcoS + /x/x^ 

(Iif2 + Kk)[Dko{T, 0 S) + co(/x/x^ 0 S + S 0 /x/x^ 
+ D{ko — Cq) vec S vec^ S 


Proof. This is obtained from Theorem 13.41 

Corollary 3.4. In particular ifY'^ ’^(/X,S 0 l£,). Then, cq = kq = 1 


E{B) = DT, + pupd' 

Cov(vecB) = (1^2 + Kif)[Zl(S 0 S) +/x/xj 0 S + S 0 /x/x^]. 


□ 

...,yD are 


□ 


)] 

□ 

, and thus 
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3.3 Method-of-moments estimators 

Returning to our notation, for which, rewrite, 0 = and /X = /X*. 

Our target is to find the method-of-moments estimators of the parameter matrices 
and /i*. First, note that the first two sample moments estimators of B are given by 

1 " 

£;(B) = - ^ B, = B = ih,), *, j = 1,..., K, 

n 

and 

Cov(vecB) = — (vecBf — vec£'(B))(vec Bf — veci?(B))^ = S. 
n 

i=l 

where S = (s^^), t,r = 1,2,..., . In addition note that for i < j and M = = 

{rriij) = and = {cij), we have 

E{b,j) = E{efBej) = efE{B)ej 

= ef {Dco'S*j^ + n* n*'^)ej = Dcgaij + rrnj, (13) 

for independent and dependent cases. 


3.3.1 Dependent case 

Note that: 


Cov{bij) = Cov(efBej) = Cov(vecef Bej) = Cov((ej 0 e^)^ vecB) 

= (ej 0 ez)"^ Cov(vecB)(ej 0 e*) 

= (e,0ed^{(Iic2+Kzc){D«o(SJ,0S;,) 

+ Co 0 Sjf + S;), 0 } 

+ D [kq — Dcq] vec vec^ 

-I- (1 — D)co [vec vec^ 

-I- vecp*/x*^ vec^ ^zc]} i^j ® s*)- 

Observing that {ej 0 ei)^Kzi' = (e^ 0 ej)'^ and that (A 0 B)(C 0 D) = (AC 0 BD), for 
i < j, ij = l,2,...,K 


Cov(6y ) = D [noaucrjj + (2ko - cl)aij] 

-I- Co [mjjau -I- muajj + 2(2 - D)mijaij]. (14) 

From (fT3l) . by replacing rriij = bij — Dcoaij in (fl^ we have 

Cov(6jj) = D{ko - 2cl)auajj + D(2 ko - {I+ 2{2 - D))cl)a'^^ 

-j- Cq \pjj^ii “t“ bii(7jj -f 2(2 E^bijCTij^ . (1^) 


Therefore equaling (TTKll to Sij = Cov{bij) we have: 


Theorem 3.6. Assume that B ~ C/PyV|^.(Zl, , Izz, D, h). Then, the method-of-moments 
estimators of and M = n* fi*^ are given by the following exact expressions. 

For i = 1,2,..., K: 




'JQ’ii “I" Sii — Qii 
TP 


( 16 ) 
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where Qf^+APsu >0,P = D{ko — 2cI) + D{2ko—{1 + 2{2 — D))cI), and Qu = 2co{3 — D)bu. 

where an has been previously found in m- 
If P — 0; thcTl (Jii — 

For i <j, i = - 1), j = 2,..A^' 


^ ^{2-DYclbl-R{T,,-s.,) - (2 - D)cob,, 

^ 

where (2 — D)‘^CQ\j — R{Tij — Sij) > 0, R = D{2 ko — (1 + 2(2 — D))cq), and 
Here an and ajj were previously computed in m- 

rriij — bij Dc^cTij ^ 


(18) 


(19) 


Denote the solution as __ ^ 

(M,S;,). 

Note that stj = Cov{bij), this is Sij are obtained from the diagonal of matrix S G 91*” . 

If R = 0, then 5^ = (s^ - Rj) / (2(2 - D)cJ>ij). 

Remark 3.2. Special attention must be payed on the constants P, Qu, R and and the 
sign of the square root, according to the selected model and the sample statistics Sy and 

bij. 


3.3.2 Independent case 

For this case, 

Cov{bij) = Cov(efBej) = Cov(vecefBej) = Cov((ej 0 e^)^ vecB) 
= (ej 0 Bi)'^ Cov(vecB)(ej 0 e^) 

= {e,(S)e,f{{lK2+KK)[DKo{'S*j,(8)-E*j,) 

+ co(p>*^ 0 + SJf 0 /X>*'^)] 

+ D{ko - Cg) vec vec^ }(ej 0 e*). 

Hence 


CoY{bij) = D [Koaiiajj + (2 ko - cl)afj] + cq {mjjau + rn^ajj + 2m^jaij). (20) 

From (fT51) . by substituting = Sij — Dcoaij in (1^01) we have 

Covi^b^j^ — D(^I‘vq 2cq^(Tu( 7jj -\~ D(2 kq 3cQ)fj^j -t- cq Ipjj^a ~t“ ^a^jj ^bijCJij^ . 

Summarising 

Theorem 3.7. Assume that yd ~ ^^(/x^,S;h), independently, for d = 1,2,..., D, 

such that 

Y = (yi|y 2 | • • • Iyd) and p = (pi|p 2 l ’ ’' IMd), 
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and let 


D 


B = YY^ = ^ydyJ. 


d=l 


Then, the method-of-moments estimators ofT,*j^ and M = ji* are given by the following 
exact expressions. 

For i = 1,2,..., K: 

~ _ Q'ii + ‘^PSii ~ Qii 


2P 

where Q% + APsa > 0, P = D{3ko — Scg), and Qu = Acobr, 

rna — Dcoctu, 

where an has been previously found in m- 
If P — 0; then On — ^nlQa. 

For i <j, i = l,...,{K - l),j = 2,...,K: 


( 21 ) 


( 22 ) 




HifTj^j ^ij') 


R 

_2 

where — R{Tij — Sij) >0,R = D{2 kq — Scg) and 


(23) 


Tij — -O(aC0 CQbjjOii ~t“ CobiiOjj . 

Here ou and ojj were previously computed in \21\) . 

niij — bj^j DcQOij, (^4) 

Denote the solution as _ 

(M,S;,). 

Note that Sij = Cov{bij), this is Sij are obtained from the diagonal of matrix S G 91^ . 

If R — 0, then o^j — ^ij ) / . 

Remark 3.3. Recall that the method-of-moments estimators are not uniquely defined. In 
addition, if instead of estimating the parameter 0, method-of-moments estimator of, say, 
g{0) is desired, it can be obtained in several ways. One way would be to first find method- 
of-moments estimator, say 6 oi 6 and then use g{9) as an estimator of g{0). Alternatively, 
we can found the moments of function g{9) and then apply the method of moments to find 
the method-of-moments estimator g{9) of g{9). Estimators using either way are termed 


meth od-of-moments estimators and may be not be the same in both cases, see (|Mood et al 


1974 Section 7.2.1, p.276). 


_ T he fol lowing result formalise the algorithm (Principal Coordinate Analysis, collected at 

Lelel ( 19931) 1 for obtain /x*, the estimated coordinates of the mean form (up to translation. 


rotation, and reflection transformations) using the method-of-moments estimator M. 

Theorem 3.8. Let M the method-of-moments estimator o/M = /x*/x*^ (for dependent or 
independent cases). LetM = ViLV^ is nonsingular part of its spectral decomposition, where 
Vi is a semiorthogonal matrix, Vi G i.e. VfVi = and L = diag(Ai,..., Au), 

with D the rank of matrix M. Then the method-of-moments estimator of fx* is 

Ji* = ViW, 

where W = diag(-yAi',..., V^d)- 
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Proof. It is follow from Remark 13.31 □ 

Theorem 3.9. Let {Ji the method-of-moments estimators of 


Then as n ^ oo 


{/j,*, (/!*, S^) in probability. 


Proof. This follows from the consistency of th e sample m oments and the continuity of the 


function (/x*,S^) in (^^(B), Cov(vecB)), see ( Raol 1973 . Section 5d.l, p. 351). 


□ 


4 Consistent estimation when T,£) is a general non-negative 
definite matrix 


Results in this section are motivated in the result obtained by Dutilleul 1 1999ll under a 
matrix multivariate Gaussian distribution via the maximum likelihood estimation. We make 
an heuristic evaluation of the useful of these results in our approach based in method-of- 
moments estimation. 

Our algorithm is based in the following modified expressions: 


1 




2=1 


Sk = 


nD 




*\T 


i=l 


Algorithm 

Initialisation: 


= 0; = SJ,; SJ, = — ^ (X^ - fsff (SJf)“ (X? - /i*); 


r = r + 1 


■^K 


1 ^ 


nD 


2=1 


= S') (x“- n ’); 


>£i or >£ 2 , 


While 

ll^r+l V'' 

W^D ^ 

Repeat: 

r = r + 1; 

’ 

\^r _ y^r-|-l. 

5 

Recompute and . 

Solutions are: 

S =*: _ V' _ V'7’ 

K — — ^D- 

Where ei and £2 define two infinitesimal positive quantities and 
norm, 1 1 IAII 2 = </tr(AA^ 


(25) 

(26) 


1 2 is the Euclidean 
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Theorem 4.1. Let (/x*, 0 ^d) the method-of-moments estimators of (SJf, /x* 0 

Then as n ^ oo 


in probability. 

Proof. This follows from Remark 13.31 


□ 


5 Estimation of the form difference 


A detailed discussion of Euclidean Distance Matrix, matri x form, form diffe rence and their 
probabilistic, geometrical, etc. properties may be found in lhelel ( 199ll [l993[l . For your con¬ 
venience, next we shall introduce some notation, although in general we adhere to standard 
notation forms. 

Consider the following square symmetric matrix, know as Euclidean Distance Matrix: 


F(X) = 


/ 0 d(l,2) 

d(2,l) 0 

Vd(i^,l) diK, 2 ) 


d{l,K-l) d{l,K)\ 
d{2, K-l) d{2, K) 


d{K,K-l) 0 


where d{i,j) denotes the Euclidean distance between landmarks i and j, in shape theory 
such matri x is termed form matrix. Among others interesting properties of form matrix, 
Leld ( 199ll i proves that F(X) is a maximal invariant under the group of transformations 


consisting of translation, rotation, and reflection. Therefor, F(X) retains all the relevant 
information about the form of an object. 

Let Xi, X2,..., X„ be n independent observation from population I and Yi , Y2 ,..., Y^ 
be m independent observation from population II. Let the mean form of population I be /x^ 
with the corresponding form matrix F (/x^) and corresponding parameters for population 
II be fjy and F (/x''^). From Leld ( I993I) we have the following definition: 

Definition 5.1. Form difference between population I and II is defined as 

FDM(/x^,/x^)=F(/x^)*F(/xY)-^, 

where * denotes the Hadamard product, 0/0 = 0 and denotes the inverse of A with 
respect to the Hadamard product, a formula for such inverse in terms of the usual product 
is given in Caro-Lopera et all ( 2012[) . 


From remark 13.31 the following theorem shows that the form difference between two 
populations can be estimated consistently when landmarks are perturbed dependently along 
each axis but independently or not correlated between the axes. 


Theorem 5.1. Let 


Kyi' 


-‘DX 


and (/x’^, S 


KY 


■‘DY 


) be the parameters for the 


two populations. If^ox — ^dy =^ 0 , then 


FDM (jl^, = F * F ^ FDM (/x^, /x^) 


in probability. 


Theorem 5.2. Let (/x^,S|fx' 

two populations. Then 


-‘DX 


and , S 


KY 


■‘DY 


) be the parameters for the 


^ ^ FDM(/x^,/x 


in probability. 
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6 Example 


The m ouse vert ebra problem was orig i nally studied in the Gaussian case bv iDrvden and Mardia 
( 1998f) (see also iMardia and DrvdenI dlQSOli'). A further analysis under elliptical models was 
implemented by Diaz-Garcia and Caro-Lopera ( 2012bll . The experiment considers the sec¬ 
ond thoracic vertebra T2 of two groups of mice: large and small. The mice are selected 
and classified according to large or small body weight; in this case, the sample consists of 
23, 23 and 30 large, small and control bones, respectively. The vertebras are digitised and 
summarised in six mathematical landmarks which are placed at points of high curvature, 
see figure [U they are symm etrically selected by measuri ng the extreme positiye and negative 
curvature of the bone. See IPrvden and Mardia ( 1998h for more details. The shape differ¬ 
ence analysis among the three groups is quite solved by a different approaches. However the 
correlation structure among landmarks requires more analysis; strong assumptions about 
those relations are usually considered because the complex exact shape distribution and a 
non existence theory for estimation for such invariant functions. 

More than an example, this landmark data is highly valuable for a correlation structure 
analysis because the symmetry of the vertebra, certainly suggest a priori a non isotropic 
model. The control group is also useful for comparisons and correctness. 


Small, large and control mouse vertebra sample 



-100 -50 0 50 too 


Figure 1: Mouse vertebra sample 


Theorems 13.61 and 13.81 can be easily implemented for a number of models. We focus on 
the main novelty (Theorem 13.61) and Kotz type model (including Gaussian) which is very 
flexible and meaningful for various values of the parameters r, s and iV, see appendix. 

First of all we illustrate Theorem 13.71 under six different models with independent land¬ 
marks. Moment-method estimates of mean shape by using the common Gaussian model is 
shown in figure [21 in this case the estimate are complete unrealistic, as we expect, given the 
assumption of independence of landmarks. However, if we consider more complex models 
based on independence, the estimation tends to be more similar to the structure suggested 
by the sample. The addressed evolution from Kotz 1 to Kotz 5 is depicted in figures [3] to |7l 
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Gaussian: S(red), L(blue), C(green) 



Figure 2: Moment method estimates under independence: Gaussian model 


Kotz 1: S(red), L(blue), C(green) 



Figure 3: Moment method estimates under independence: Kotz 1 model 


An heuristic behavior is noted, the lack of dependence in the Gaussian model, and its 
unrealistic moment method estimates, it seems to be improved by considering a more robust 
Kotz type model even with landmark dependence. The literature has studied this artificial 
data by the independent Gaussian case, so, given that no expert have set this assumption 
we can get further into more robust analysis and advance in some selection criteria, but if we 
have an experiment modeled by the independent Gaussian case according to the opinion of 
an expert in the field, we must follow that law and the further analysis, based on landmark 
dependence and elliptical families, that we provide next, cannot be implemented in such 
cases. 

In the artificial mice data, we now can focus on the dependent case and the moment 
method estimators of Theorem l3.61 given that the Gaussian case is out of any consideration, 
then we have to study, for example, other Kotz models. In order to illustrate the important 
effect of landmark dependence we consider the simplest Kotz model after Gaussian, when 
iV = 2, r = 1/2 and s = 1, which is referred as Kotz 1 model, and we compare the perfor- 
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Kotz 2: S(red), L(blue), C(green) 



Figure 4: Moment method estimates under independence: Kotz 2 model 


Kotz 3: S(red), L(blue), C(green) 



-0.4 -0.2 0.0 0.2 0.4 


Figure 5: Moment method estimates under independence: Kotz 3 model 


mance of Theorem 13.61 with another mean shape estimations. Table [2] provides comparisons 
among mean shape estimates of the small group, the y include me an shape by moments of 
The orem 13.61 the mea n shape by Frechet method fsee lKend ( 1992r) '). and Bookstein method 


Bookstein ( IQSdlB : certainly the estimations are truly similar. Note also that Kotz 1 


(_see 

law with independent landmarks provided a bad moment method estimator of mean shape, 
but the same model under the expected and realistic dependence revels similarity with more 
complex mean shape estimators derived by standard shape theories, see figures [3] and [3 
respectively. 

The exact formula for the moments estimation ITheorem 13.61) also agrees with the pre¬ 
vious conclusions in literature about strong difference in Gaussian mean shape between the 
small (S) and large (L) groups. Figure [3 also shows the mean shape estimation of the con¬ 
trol (C) group. As we expect, the control group must tend to show strong symmetry among 
landmarks, by ’’averaging” in some sense the small and large estimates. 

Different types of Kotz distribution have also modeled the sample, they correspond to 
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Kotz 4: S(red), L(blue), C(green) 



- 0.4 - 0.2 0.0 0.2 0.4 


Figure 6: Moment method estimates under independence: Kotz 4 model 


Kotz 5: S(red), L(blue), C(green) 



Figure 7: Moment method estimates under independence: Kotz 5 model 


Table 2: Estimation for the mean shape for the small group by Theorem 13.61 (Kotz 1), 
Frechet (F), and Bookstein (B). 


Th. 6, 

^oTs 

0.5 

0.084507028 

0.014836162 

- 0.073397569 

- 0.005026754 


Th. 6, M2 
0 
0 

0.3301634 

0.6957339 

0.3394693 

- 0.2184060 


B. Ml 
- 0.5 
0.5 

0.08469746 

0.01215768 

- 0.06874750 

- 0.02502185 


B. M2 
0 
0 

0.2933430 

0.5613175 

0.2991278 

- 0.3041418 


F. Ml 
- 0.5 
0.5 

0.08490820 

0.01245608 

- 0.06869796 

- 0.02512807 


F- M2 
0 
0 

0.2924684 

0.5589496 

0.2982314 

- 0.3044915 


the denoted models Kotz 1, Kotz 2, Kotz 3, Kotz 4 and Kotz 5, with parameters N = 
2,s = l,r = 1/2; iV = 3,s = l,r = 1/2; N = 2,s = 2,r = 1/2; iV = 2, s = 3, r = 1/2 
and N = 20,s = 20,r = 1/2, respectively. Technical details about the generalised singular 
Pseudo-Wishart distributions and particular Kotz Pseudo-Wishart distributions referred in 
this example, can be seen in the appendix. The corresponding mean shapes estimates were 
computed, but for reasons of space, we only show the results of the Kotz 5 model (suggested 
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Kotz 1 (Th. 6): S(red), L(blue), C(green) 



Figure 8: Moment method estimates under dependence: Kotz 1 model 


by the preceding independent results and certain selection criteria that we will propose later) 
see figure [H 


Kotz 5 mean shape by Theorem 6 



Figure 9: Mean shape estimates for large, small and control groups under the Kotz 5 model 


Now we apply the algorithm for a consistent estimation when is a general non¬ 
negative definite matrix, under the Kotz 5 law. For the routine propose in Section|T]we have 
fixed El = £2 = 0.000005 in the three groups small, large and control, then we found that 
the number of iteration to reach the addressed tolerance is 57, 53 and 61, respectively. 

For the small group the estimated covariance matrices are given next (here the asso¬ 
ciated correlation matrix p of S is provided, for the sake of interpretation, recall that 
p= (diag(S))~^ S (diag(S))”5): 


/ 1.0000000 

-0.88123087 

-0.45286210 

-0.0947470 

0.3167573 

0.1049686 

\ 

-0.8812309 

1.00000000 

-0.01931031 

-0.3859582 

-0.7246992 

0.3741399 


-0.4528621 

-0.01931031 

1.00000000 

0.9267571 

0.6990041 

-0.9323209 


-0.0947470 

-0.38595825 

0.92675709 

1.0000000 

0.9113738 

-0.9979133 


0.3167573 

-0.72469917 

0.69900414 

0.9113738 

1.0000000 

-0.9087033 


\ 0.1049686 

0.37413987 

-0.93232089 

-0.9979133 

-0.9087033 

1.0000000 

/ 
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and 


1.00000000 -0.1305434 

-0.1305434 1.00000000 


Pd — 


( 


For the large group the estimated correlation matrices are: 


/ 

1.00000000 

-0.65790585 

-0.49909037 

-0.05610956 

-0.07744806 

0.42376772 



-0.65790585 

1.00000000 

0.06499537 

-0.44441269 

-0.29572585 

0.03327717 



-0.49909037 

0.06499537 

1.00000000 

0.42647926 

0.65775361 

-0.76574318 

Pk = 


-0.05610956 

-0.44441269 

0.42647926 

1.00000000 

0.65493138 

-0.73668655 



-0.07744806 

-0.29572585 

0.65775361 

0.65493138 

1.00000000 

-0.79064351 


1 

0.42376772 

0.03327717 

-0.76574318 

-0.73668655 

-0.79064351 

1.00000000 

and 












^ 1.00000000 

-0.2080039 \ 






Pd — 1 

-0.2080039 

1.00000000 i 




Meanwhile in the control group the estimated correlation matrices are: 


I 

Pk = 

\ 

and 


1.00000000 

-0.61798019 

-0.61373047 

-0.39687003 

0.05198526 

0.50362856 


-0.6179802 

1.0000000 

-0.1036402 

-0.4037052 

-0.7081139 

0.2660900 


-0.6137305 

-0.1036402 

1.0000000 

0.7604036 

0.5048809 

-0.8386367 


-0.3968700 

-0.4037052 

0.7604036 

1.0000000 

0.6481050 

-0.9587295 


0.05198526 

-0.70811391 

0.50488092 

0.64810504 

1.00000000 

-0.64271257 


0.5036286 \ 

0.2660900 
-0.8386367 
-0.9587295 
-0.6427126 
1.0000000 / 


f 1.00000000 0.1048453 \ 

1^ 0.1048453 1.00000000 ) ' 


The three groups reveal almost null correlation among axes, but some important correla¬ 
tion among landmarks, as we expect from the pseudo-symmetry of the bones. The estimates 
in the small and large groups detects the main landmarks responsible for the mean shape 
difference, meanwhile in the control case the estimates tends to follow the main contribution 
of large or small differentiating landmarks as we expect. 


In a similar way we have run the routines with the same tolerance ei = £2 = 0.000005 
for the models Kotz 1, to Kotz 4; they reached the stability between 50 to 70 iterations 
in the three groups, and similar conclusions about the almost null correlation among axes 
and strong correlation among landmarks were found in the models. We will not show the 
estimates of each Kotz type, but have provided the results for the model Kotz 5 type, for 
reasons that we will explain later when the ’’best” model is selected under certain criteria; 
it was also suggested by the independent case analysis. 

For a selection model criteria, the control group plays a fundamental role, in this case we 
just need to look for the law which obtains the minimum coefficient of variation when the 
small and large groups are compared with the control one; the analysis also must consider 
the distance between small and the large group relative to the mean wi th controls. W e appl y 
non-Euclidian distance between covariance matrix, a technique due to Drvden et al.l ( 2009ll . 
The method is approjniate for meaningful correlation matrices, in this case it is performed 
only for , because certainly ratifies in all the models that no correlation among axes 
is observed. In Tables [3] and 01 Kl,..., K5, s, I, c, stand for Kotz 1,..., Kotz 5, small, large 
and control, respectively. 


Tables [3] and [4] shows all the pairwise covariance distances, in particular, the percentage 
variation coefficient is presented in parenthesis. We are searching for models which reflect 
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Table 3: Selection model criteria. 



Kll 

Klc 

K2s 

K21 

K2c 

K3s 

K31 

K3c 

Kls 

12.9 

8.7(37) 

12.8 

15.9 

14.1 

11.8 

15.7 

13.7 

Kll 


5.1(37) 

10.6 

6.1 

8.4 

10.8 

5.1 

8.5 

Klc 



9.6 

9.1 

8.9 

9.2 

8.6 

8.7 

K2s 




11.0 

11.0(14) 

6.0 

11.1 

10.7 

K21 





9.0(14) 

12.0 

1.5 

9.1 

K2c 






11.6 

8.8 

0.9 

K3s 

K31 







12.1 

11.2(15) 

9.0(15) 


Table 4: Selection model criteria. 



K4s 

K41 

K4c 

K5s 

K51 

K5c 

Kls 

11.1 

15.0 

13.1 

11.1 

13.9 

12.1 

Kll 

11.3 

3.9 

8.6 

12.2 

2.4 

2.8 

Klc 

9.2 

7.7 

8.4 

9.7 

6.4 

4.8 

K2s 

8.9 

10.8 

10.3 

10.9 

10.1 

9.7 

K21 

13.3 

3.0 

9.5 

14.7 

4.4 

7.0 

K2c 

12.4 

8.5 

2.1 

13.4 

8.2 

8.4 

K3s 

6.4 

11.6 

10.7 

9.2 

10.7 

9.8 

K31 

13.2 

1.6 

9.5 

14.6 

3.3 

6.2 

K3c 

11.9 

8.7 

1.2 

12.9 

8.4 

8.3 

K4s 


12.6 

11.4(16) 

6.8 

11.6 

10.6 

K41 



9.1(16) 

13.9 

1.8 

5.1 

K4c 




12.3 

8.7 

8.3 

K5s 





12.8 

11.6(74) 

K51 






3.6(74) 


the role of the control group and separated the classes properly, the analysis must be com¬ 
plemented with mean shape estimates and a third criteria involving how the last ones are 
far from another accepted estimates, the Frechet mean shape for example. The a ddressed 
mean shape distance can be achieved by a number of approaches, see for example Kendalll 
(Il984ll . 


Kotz 3 and 4 laws behave well with percentage variation coefficient, but the corresponding 
distance with the control group and the sample is to far to be realistic, specially with the 
small group. If we find the so called Riemannian distance among the moment method 
estimates and Frechet and Bookstein mean shape, we obtain the results of tabled 


Table 5: Model selection criteria. 



K2 

K3 

K4 

K5 

F. 

B. 

K1 

0.274 

0.236 

0.211 

0.180 

0.113 

0.112 

K2 


0.082 

0.128 

0.153 

0.189 

0.191 

K3 



0.048 

0.078 

0.131 

0.133 

K4 




0.035 

0.099 

0.100 

K5 





0.067 

0.068 

F 






0.002 


The mean shape estimate based on a Kotz 5 model is very near to the estimates computed 
by Frechet and Bookstein (which are significantly similar), it also reflects good difference 
between the small and large; similar findings for the large and small group were computed, 
then collecting the results, we can propose Kotz type 5 model as a suitable law for modeling 
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this particular example. Note that the selection agrees with the conclusion proposed in the 
independent case. 

It is important to note, that mathematical or statistical selection is just a suggestion for 
an experiment lacking of any prior assumption of the supporting distribution provided by an 
expert. In our case, literature shows no ex pert assumption about normality, in fact, this data 
full studied in iDrvden and Mardial (1998) and the references therein, was traditionally set in 
the Gaussian theory in order to simplify computations and/or the use of the classes of exact 
distributions were not available at that time. However, if an experiment was sufficiently 
studied by an expert which the Gaussian model is truly normal, then the above selection 
of models, are out of significance; and given that moments-method estimates does not work 
under Gaussianity, as we have shown in this example, then the results presented here cannot 
be applied properly. 

Now, at this stage, the conclusion about Kotz 5 model ratifies that non-Gaussian models 
explain better t he three samples (an elliptical isotropic a pproach also verified this conclusion, 
see for example Dfaz-Garcfa and Caro-Lopera ( 2012bll ). 

Once the model is selected, we are interested in application of Section^ about estimation 
of mean form difference. In fact, we can go further by considering hypothesis testing for 
equality of the associated Euclidean Distance Matrices of t wo po pulations. 

The methodology can be found in lLele and Richtsmeierl ( I99lll and the references therein. 
We are interested in testing Hq : F(/x^) = cF(/i'^), for some c > 0, where and are 
the population mean shape. Based on a sample of objects X’s and Y’s, with corresponding 
estimated mean shapes and obtained with the exact formula given in Theorem 


we derive the form difference matrix FDM 




This last matrix can be used for 


defining a number of suitable statistics for testing Hq, however. iLele and Richtsmeien (1199111 
recommend the following: 


T = max FDMij j min FDM, 






i] I 


where FDMij is the element of matrix FDM. Note that if F[q is tru e T is close to one. 
Moreo ver, T satisfies the desirable property of invariance under scaling, see lLele and Richtsmeiei 
(|l9911) for more details. 

The null distribution is difficult to obtain even in the simplest case of Gaussian, so we 
can obtain an empirical nul l distribution by using the well known bootstrap procedure, see 
Lele and Richtsmeien (1199111 and the references therein. For similar samples of the current 
example, the above referred authors recommend a bootstrap of size 100. 

Once the empirical distribution is obtain, a p-value, based on the upper tail of the 
observed statistics, rejects FIq for small values near to 0.1. 

Table [S] reports such tests for the Gaussian and Kotz 3 type models and the three 
pair comparisons of interest. We note that the usual Gaussian case, under the expected 
dependence condition of Theorem 13.61 malfunction and cannot detect the role of the control 
test, given a wrong conclusion. The selected model by covariance distances, separates as 
we expect the control group and gives a suitable p-value of certain difference, but it is 
not sufficient enough for concluding shape difference. This open an interesting dis cussion 
about the method based on coordinate free approach of Lele and Richtsmeier ( 199lll . given 
that the quotient pairwise-element in definition of the matrix form difference is neglecting 
the whole matrix structure. Improving this aspect deserves a further work by defining 
a more robust matrix FDM based on usual products than the very restrictive Hadamard 
product. Moreover, finding the corresponding exact distribution of FDM(X, Y) can provide 
a promising null distribution which can model hypothesis testing efficiently. 
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Table 6: p-values of testing the equality of mean shape under different models and pairs of 
populations . 



Small-Large 

Small-Gontrol 

Large-Gontrol 

Gaussian 

0.00 

0.00 

0.00 

Kotz 5 

0.12 

0.51 

0.74 


7 Conclusions 


1. First, by replacing the Gaussian model for the elliptical model, an infinite range of 
possibilities in making an assumption of a model is opened, allowing to model a wide 
range of real situations, more or less heavy tails and more or less kurtosis than the 
Gaussian model. 


2. Under this family of elliptical models is possible consistently estimate all parameters. 

3. Notably, all these estimators are extremely easy to calculate. 

4. Alternatively to the hypothesis assumed in subsections 2.3 and 2.4, an interesting 
alternative to investigate is: Assume that the joint distribution of Ei, E 2 ,..., E„ is 


E = (El, E 2 ,..., E„) ~ SKxnoiO, S 


K ' 


-‘D 


’ In 5 ? 


where Gov(vecE^) = S^0S p0ln . Then proceeding in a s imilar way and generalising 
to no central case the r esults in ( Fang and Zhang . 1990l Eq. 3.4.14, p. 109) and 
( Gupta and Varg£] . ll993L Theorem 5.1.6, ). 170), we obtain that the joint distribution 


of Bi,B 2 ,...,B„ is 


B = (Bi,B2,...,B„)^01PWK.n , 

where fl = (S)^)” Noting that in this case it is implicitly assuming that 

the sample Xi,..., X„ is dependent. 

5. Alternatively, and recalling that the method-of-moments estimators are not uniquely 
defined (see Remark 1531) the method-of-moments estimator of Sd can be obtained 
from first two moments of B too. 


Acknowledgments 

This article was written under the existing research agreement between the hrst author and 
the Universidad Autonoma Agraria Antonio Narro, Saltillo, Mexico. F. Garo was supported 
by a research project of University of Medellin, Golombia. 

References 

Arnold, S. F. (1981). Theory of linear models and multivariate analysis, John Wiley & 
Sons, Inc., New York. 

Bookstein, F. L. (1986). Size and shape spaces for landmark data in two dimensions (with 
discussion), Stat. Sci. 1 181-242. 


29 




















Caro-Lopera, F. C. and DfAZ-GARcfA, J. A. (2012). Matrix Kummer-Pearson VII 
relation and polynomial pearson VII configuration density, J. Iranian Stat. Soc. 11(2) 
217-230. 

Caro-Lopera, F. J., Diaz-Garcia, J. A., and Gonzalez-Farias, G. (2009). Noncen¬ 
tral elliptical configuration density, J. Multivariate Anal. 101(1) 32-43. 

Caro-Lopera, F. J., Diaz-Garci'a, J. A., and Gonzalez-Fari'as, G. (2014). Inference 
in affine shape theory under elliptical models, J. Korean Stat. Soc. 43(1) 67-77. 

Caro-Lopera, F. J., Leiva, V. and Balakrishnan. N. (2012). Connection between the 
Hadamard and matrix products with an application to matrix-variate Birnbaum-Saunders 
distributions, J. Multivariate Anal. 104 126-139. 

DfAZ-GARcfA, J. A., AND Caro-Lopera, F. J. (2012a). Generalised shape theory via SV 
decomposition I, Metrika 75(4) 541-565. 

Diaz-Garcia, j. a., and Caro-Lopera, F. J. (2012b). Statistical theory of shape under 
elliptical models and singular value decompositions, J. Multivariate Anal. 103(1) 77-92. 

Diaz-Garcia, J. A., and Caro-Lopera, F. J. (2013). Generalised shape theory via 
pseudo Wishart distribution, Sankhya A 75(2) 253-276. 

Diaz-Garcia, J. A., and Caro-Lopera, F. J. (2014). Statistical theory of shape under 
elliptical models via QR decomposition. Statistics 48(2) 456-472. 

Diaz-Garcia, J. A. (1994). Contributions to the theory of Wishart and multivariate ellip¬ 
tical distributions, Ph.D. disertation, Universidad de Granada, Spain, (in Spanish). 

DIaz-GarcIa, j. a., and R. Gutirrez Jimez, R. (1996). Matrix differential calculus 
and moments of a random matrix elliptical, Serie Coleccin ’’Estadstica Multivariable y 
Procesos Estocsticos”. Universidad de Granada, Espaa, (in Spanish). 

Di'az-Garci'a, j. a., and Gonzalez-FarIas, G. (2005). Singular Random Matrix de¬ 
compositions: Distributions, J. Multivariate Anal. 94(1) 109-122. 

Diaz-Garcia, J. A., Gutierrez, R. J. and Ramos-Quiroga, R. (2003). Size-and-shape 
cone, shape disk and configuration densities for the elliptical models, Brazilian J. Prob. 
Stall. 17 135-146. 

DfAZ-GARCfA, J. A., AND GuTiERREZ-JAIMEZ, R. (2006). Wishart and Pseudo-Wishart 
distributions under elliptical laws and related distributions in the shape theory context, 
J. Statist. Plann. Inference 136(12) 4176-4193. 

Dryden, I. L. AND AND Mardia, K. V. (1998). Statistical shape analysis, John Wiley and 
Sons, Chichester. 

Dryden, I. L., Koloydenko, A., and Zhou D. (2009). Non-Euclidean statistics for 
covariance matrices, with applications to diffusion tensor imaging. An. App. Statist. 3 
1102-1123. 

Dutilleul, P. (1999). The mle algorithm for the matrix normal distribution, J. Statist. 
Comput. Simul. 64 105-123. 

Fang, K. T. and Zhang, Y. T. (1990). Generalized multivariate analysis, Science Press, 
Beijing, Springer-Verlang. 


30 


Fang, K. T., Kotz, S., and Ng, K. W. (1990). Symmetric multivariate and related 
distributions, Chapman Hall, London. 

Goodall, C. R. (1991). Procustes methods in the statistical analysis of shape (with dis¬ 
cussion), J. Royal Stat. Soc. B 53 285-339. 

Goodall, C. R., and Mardia, K. V. (1993). Multivariate aspects of shape theory, Ann. 
Statist. 21 848-866. 

Gupta, A. K., and Varga, T. (1993). Elliptically contoured models in statistics, Kluwer 
Academic Publishers, Dordrecht 

Kent D. G. (1984). Shape manifolds, Procrustean metrics and complex projective spaces, 
Bulletin of the London Mathematical Society, 16 81-121. 

Kent J. T. (1992). New directions in shape analysis. In Mardia, K. V., editor. The Art of 
Statistical Science, 115-127. Wiley, Chichester. 

Khatri, C. G. (1968). Some results for the singular normal multivariate regression nodels, 
Sankhyd A 30 267-280. 

Koev, P., and Edelman, a. (2006). The efficient evaluation of the hypergeometric func¬ 
tion of a matrix argument. Math. Comp. 75 833-846. 

Le, H. L., and Kendall, D. G. (1993). The Riemannian structure of Euclidean spaces: 
a novel environment for statistics, Ann. Statist. 21 1225-1271. 

Lele, S. (1991). Some Comments on Coordinate Free and Scale Invariant Method in Mor¬ 
phometries, Am. J. Phys. Anthropol. 85 407-418. 

Lele, S. (1993). Euclidean distance matrix analysis (EDMA): Estimation of mean form and 
mean form difference, Math. Geol. 25(5) 573-602. 

Lele, S., and Rightsmeier, J. (1991). Euclidean distance matrix analysis: A coordi¬ 
nate free approach for comparing biological shapes using landmark data, Amer. J. Phys. 
Anthropol. 86 415-427. 

Lele, S., and Rightsmeier, J. (1990). Statistical models in morphometries: Are they 
realistic? Syst. Zool. 39(1) 60-69,. 

Magnus, J. R., and Neudecker, H. (1979). The commutation matrix: Some properties 
and applications, Ann. Statist. 7(2) 381-394. 

Mardia, K. V., and Dryden, 1. L. (1989) The Statistical Analysis of Shape Data, 
Biometrika 76 271-281. 

Mood, A. M., Graybill, F. A., and Does, D. G. (1974). Introduction to the theory 
statistics. Third edition, McGraw-Hill Series in Probability and Statistics New York. 

Muirhead, R. j. (1982). Aspects of multivariate statistical theory, Wiley Series in Proba¬ 
bility and Mathematical Statistics. John Wiley & Sons, Inc., New York. 

Nadarajah, S. (2003). The Kotz-type distribution with applications. Statistics 37(4) 341- 
358. 

Rao, C. R. (1973). Linear Statistical Inference and Its Applications, John Wiley & Sons, 
New York, 1973. 


31 


Richtsmeier, J. T., Deleon, V. B., and Lele, S. R. (2002). The Promise of Geometric 
Morphometries, Yearbook of Phis. Anthropol. 45 (2002) 63-91. 

Walker, J. (2001). Ability of geometric morphometric methods to estimate a known co- 
variance matrix, Syst. Biol. 49 686-696. 


A Particular generalised Pseudo-Wishart singular dis¬ 
tributions 


The fo llo wing result is a particular case of the gener al result in iDiaz-Garcia and Gonzalez-Farias 


( 2005 ) or Diaz-Garcia and Gutierrez- Jaime3 ( 20061) , when © is non singular and the notation 
of this paper is assumed. 

Theorem A.l (Generalised singular Pseudo-Wishart distributions). Assume that Y 
^KxD ^S G 0, h), where h admits a power series expansion 


h{v + a) =Y^ 


h^'^\a)i 


t=o 


in 5R. Let, also, q = min(A" — 1,-D); then the density o/B = Y© ^Y^ is given by 

^gD/2|L|(D-A-l)/2 “ h(2‘)(tr(S-B + U)) ^^(f^S-B) 

■2^2^-n-mTR-(^B) 






D/2 


t—0 K 




(27) 


where B = WiLW^, is the nonsingular espectral decomposition ofB withWi a semiorthog- 
onal matrix, i.e. W^Wi = Ip , and L = diag(/i,..., G); = Y, (dB) is 

Hausdorff measure is defined in IPiaz- Garcia and Gonzalez-Farias \ . \200A . Section 5); \i, 
i = 1,..., (AT — 1), are nonull eigenvalues of S, and where Ck,{A.) are the zonal polyno¬ 
mials of A corresponding to the partition k = {ti,... ,ta) oft, with = 

11^=1 (® ~ 0 ~ l)/2)tj 7 (a)t = a(a -I- 1) ■ ■ ■ (a -h t — 1), being the generalized hypergeometric 
co effficients an d TAa ) = 0^=1 r(a ~ {j ~ l)/2) is the multivariate gamma function. 


see 


Muirhead \ 


Corollary A.l (Singular Pseudo-Wishart Gaussian distribution). Let us suppose thatY ~ 
Mkxd (/^’ ^ ® ®)' “ min(A — 1, D); then the density o/B = Y©~^Y^ is given 

by 

^ 1 /I i„„ 

(28) 


= Cetr - -(S^B - n) b^i - D; -flS"B (dB), 


1 1 


with 


C = 


^r,(,_(K_l))/2|L|(D-A-l)/2 


CK-1 


2^(^-i)/2r,[iD] Y[ X 


D/2 


. i=l 


where oA'i(-) is a hypergeometric function with a matrix argument, see ItMuirhead \1982 . p. 
258). 
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B Singular Pseudo-Wishart Kotz distribution. 

Firs recall that the K x D random matrix X is said to have a singular matrix multivariate 
symmetric Kotz type distribution with parameters N,r, s € iR, : K x : K x K, of 

rank K — 1, & : D x D with r > 0, s > 0, 2N + {K — 1)D > 2, S > 0, and © > 0 if its 
density is 


g^(2N+{K-l)D-2)/2sY _ 1)D/2) 


/K-1 


^(k-i)d/2y [(2iv +{K- l)D - 2)/2s\ Af |©|(x-i )/2 


X [tr©~i(Y-/^)^S'(Y-/x)]^ ^ exp {-rtr"©-\Y-/x)'^S”(Y-//)}. 

When T = s = 1, and i? = 1/2 we get the singular matrix variate gaussian distribution. 

Note that particular singular Pseudo-Wishart distributions just depend on the general 
derivative of the elliptical generator function; it seems a trivial fact, but the gen¬ 

eral fo rmulae involves cumbersome expressions indexed by partitions, see lCaro-Lopera et al.l 
( 2009ll . In the case of Kotz type distribution they derived the following expressions. 

When s = 1, the Kotz type models and their general derivative simplify substantially. 
Thus, the following expressions applies for Gaussian, and the so called Kotz 1, Kotz 2, with 
parameters = 1, s = 1, r = 1/2; iV = 2, s = 1, r = 1/2; iV = 3, s = 1, r = 1/2; respectively. 
The generator model is given by 

( 2 ^) Tr(^-i)D/2YlN-l + {K-l)D/2f 
And, the corresponding /c-th derivative of h, follows from 

N-l r 1 

■^y exp[-r2/]. 


which is given by 


{-r)^y^ ^ exp[-r?/] 11-b ^ 

V—1 ^ ^ 


'v — 1 


n(iv-i-*) 


li^O 


{-ryY 


where k = 2t. 

For the remaining models of the example, the so termed Kotz 3, Kotz 4 and Kotz 5, 
with parameters N = 2, s = 2, r = 1/2; N = 2, s = 3, r = 1/2 and N = 20, s = 20, r = 1/2, 
respectively, the generator function is given by: 

^^i2N+iK-l)D-2)/2s r[(K - l)D/2] 

(y) Tr{K-i)D/2Y[(^2N+{K-l)D-2)/2s]^ ^ ^ 


meanwhile the required fc-th derivative of h, follows from exp (—ry®), which is given by 


T-l -Ry^ 

y e 


E 

KSPk 




k\{-R) 




wuvm- 


k 

+ E 

m=l 



m— 1 

n(T-i-i) 


. i=0 
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X E 

KGPk-, 


{k — m)!(—7?)^*=i 


n';zr\s-j)^"=r+. 


E k — 7TL f ' \ 

i = l {s-%)vi-r. 




where denotes the summation over all the partitions 

K = (fc - ..., 

of k, with J 2 i=i i-e. k is a partition of k consisting of vi ones, V 2 twos, V 3 threes, 

etc. It is important to quote that all the singular Pseudo-Wishart distributions associ¬ 
ated with the above Kotz type kernels can be co mputed by some modifications of the al- 
gorithms provided by Koev and EdelmaiTI ( 200fi[l for the Gaussian case, see for example 
Diaz-Garcia and Caro-Loperal ( 2013ll and similar works of the authors on shape theory. 
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